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The eigenstates of many-body localized (MBL) Hamiltonians exhibit low entanglement. We adapt 
the highly successful density-matrix renormalization group method, which is usually used to find 
modestly entangled ground states of local Hamiltonians, to find individual highly excited eigen¬ 
states of MBL Hamiltonians. The adaptation builds on the distinctive spatial structure of such 
eigenstates. We benchmark our method against the well studied random field Heisenberg model in 
one dimension. At moderate to large disorder, the method successfully obtains excited eigenstates 
with high accuracy, thereby enabling a study of MBL systems at much larger system sizes than 
those accessible to exact-diagonalization methods. 


Introduction: Many-body localization generalizes An¬ 
derson localization to interacting systems and entails dis¬ 
order induced breakdown of ergodicity and thermaliza- 
tion. Its existence was only recently settled, following 
precursors [1-3], via a series of perturbative arguments 
[4, 5] and numerical studies [6-10]. An intense effort has 
followed revealing an extremely rich set of properties ex¬ 
hibited by MBL systems; see e.g. the review [11] and 
references therein. This work has revealed the central¬ 
ity of many-body eigenstates to understanding a regime 
where quantum statistical mechanics simply does not ap¬ 
ply. The phase transition between the ergodic and lo¬ 
calized regimes]?, 10, 12-16] involves a singular change 
in the entanglement entropy of eigenstates from volume 
law to area law. Moreover, further transitions in the lo¬ 
calized regime can involve the development of eigenstate 
order [16-21], whereby individual highly-excited eigen¬ 
states can display patterns of order (both symmetry¬ 
breaking and topological) that may even be forbidden in 
an equilibrium setting. Such eigenstate phase transitions 
can take place at T = oo and even while usual thermody¬ 
namic functions remain non-singular. As such transitions 
are completely invisible to the traditional tools for study¬ 
ing finite energy-density phases, they necessitate a study 
of individual highly excited eigenstates. 

Since typical MBL eigenstates have only local, area 
law, entanglement [7, 22]—although deviations from the 
area law due to rare many-body resonances and Grif¬ 
fiths effects are a complication to bear in mind—the 
well known connection between area laws and matrix- 
product state (MPS)/tensor-network representations of 
many-body states [23-26] implies that they can be effi¬ 
ciently described, even at large L. Indeed, Pekker and 
Clark [27] have examined the unitary operators that ex¬ 
actly diagonalize fully MBL (fMBL) systems and shown 
that they can be represented efficiently—in contrast to 
delocalized systems which require the full many-body 
Hilbert space for their specification. Parallel work [28] 
argued for the existence of a single “spectral tensor net¬ 
work” that efficiently represents the entire eigenspectrum 


of fMBL systems. Recently the present authors and 
J. 1. Cirac developed an efficient variational algorithm 
[29] to actually find an approximate, compact representa¬ 
tion of the diagonalizing unitary for fMBL systems—and 
hence obtain all the eigenstates. This algorithm captures 
the gross features of the spectrum very well, but does not 
target individual eigenstates to high accuracy. Here we 
describe an alternative, complementary, procedure that 
can be used obtain specific excited MBL eigenstates to 
high accuracy for large system sizes. 

Our approach is directly inspired by the density ma¬ 
trix renormalization group (DMRG) [30, 31] which has 
been used to great effect to obtain modestly entangled 
ground states in low-dimensional systems. In the MPS 
formalism, the DMRG algorithm variationally optimizes 
the MPS to minimize the ground state energy of a given 
Hamiltonian H. Naively, we could modify this algorithm 
for MBL systems by targeting the eigenstate with energy 
closest to a specified excitation energy. However this is, 
in its simplest form, problematic due to the extremely 
small— 0{e ~^)—generic many-body level spacings as we 
will explicitly show below. 

We show that this problem can be overcome by making 
use of a dehning characteristic of MBL phases, namely 
the existence of an emergent set of L commuting Z 2 - 
valued local integrals of motion (often called “1-bits”) [32- 
35]. Importantly, neighboring eigenstates with respect to 
the energy differ extensively in their spatial properties— 
we must typically flip 0{L) 1-bits to go between them— 
while there is a (soft) gap to excitations with finite num¬ 
bers of 1-bit flips. This leads to a natural algorithm in 
which we select excited eigenstates based on their over¬ 
lap with particular, localized spatial patterns instead of 
their proximity to particular energies. By this overlap 
metric, “nearby” states differ by a few [0(1)] flips of lo¬ 
cal “1-bits. But such states are typically far separated in 
energy and thus the danger of mixing in eigenstates with 
exponentially small energy splittings is minimized. 

We start with a brief review of the ground state DMRG 
method before describing our modified DMRG-X proce- 
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dure. We then apply the method to the random field 
Heisenberg chain and evaluate our results using vari¬ 
ous metrics like energy variances and overlaps with ex¬ 
act eigenstates. For strong enough disorder, we obtain 
eigenstates with machine-precision variance and find a 
rapid convergence of variances with bond dimension. Fi¬ 
nally, we use our eigenstates to efficiently compute local- 
expectation values and demonstrate the failure of the 
eigenstate thermalization hypothesis (ETH) [36-38] in 
the MBL phase. We note unpublished work [39, 40] that 
also generalizes DMRG to highly excited states using a 
more complex energy based targeting approach. 
DMRG-X Method: The proposed method is a refor¬ 
mulation of the standard DMRG algorithm [30, 31] to 
find highly-excited states of MBL systems. For a one¬ 
dimensional system of L sites, a general quantum state 
jT) can be written in the following MPS form: 

14/)= ^ (1) 

31,■■■.3L 

Here, is a x „ x Xn+i matrix and \jn) with = 

1,..., d is a basis of local states at site n (for a spin 1/2 
system, d = 2). Each matrix product ]/[j HWli 
produces a complex number which is the amplitude of 
Jd/) on the basis state \ji ■ ■ ■ Jl)- The key insight behind 
the success of DMRG is that ground states of one di¬ 
mensional systems are efficiently approximated by MPS 
[26]. Starting from an initial random MPS, the DMRG 
algorithm sweeps through the system and iteratively op¬ 
timizes the matrices Ht"!-/" by locally minimizing the en¬ 
ergy with respect to a given Hamiltonian H [41]. For 
the commonly used two-site update which simultaneously 
updates the matrices and i3["+ib"+ij an effective 

Hamiltonian % is constructed by projecting H to a mixed 
XnXn+ 2 d'^ dimensional basis . Here, the local basis states 

represent the two updated sites, and the eigen¬ 
states of the reduced density matrix lxn)L and |xn-i- 2 )fl 
compactly represent the environment to the left and right 
of the updated sites. The ground state of % is found and 
the matrices on sites n, n-l-1 are updated. The procedure 
is then repeated for all sites until convergence is achieved. 

The DMRG-X method for finding excited eigenstates 
proceeds similarly to the standard DMRG algorithm in 
that we iteratively optimize an MPS. The key difference 
is that the algorithm does not attempt convergence in 
the energy of Td but instead in the local spatial structure 
of the eigenstate. We start by initializing the algorithm 
with a product state that has a finite overlap with some 
1-bit state, e.g., for the random-field Heisenberg model 
discussed below, we choose random states in the basis 
of the form = | tiit it)- We start our DMRG- 
X algorithm with the following local two-site update : 
(i) Gonstruct the effective Hamiltonian Td. (ii) Pick the 
eigenstate of Td that has maximum overlap with the cur¬ 
rent MPS. (hi) Update the tensors and 

To produce the data below, we use a full diagonalization 
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FIG. 1. Comparison between eigenenergies obtained nsing 
exact diagonalization (blue) and variational DMRG-X (red) 
for a system of size L = 12, disorder strength IF = 8 and bond 
dimension x = 16. The successive panels which zoom into 
the shaded regions of the spectrum show that all individual 
eigenenergies are obtained extremely accurately. The bottom 
right panel shows the exact and variational amplitudes for a 
particular eigenstate with Mott-resonances, showing that the 
method successively captures resonant states. 

of Td which scales as x^ with x being the mean bond di¬ 
mension. Alternatively, it is also possible to find a small 
set of k eigenstates of Td near the energy of the current 
MPS and then pick the eigenstate with largest overlap. 
This yields an algorithm that scales approximately as x^i 
but an optimal k has to be found for each case. 

This DMRG-X prescription ensures that no individ¬ 
ual update step of the MPS matrices results in a large 
spatial reorganization, which is appropriate for a local¬ 
ized phase. By contrast, if we pick excited eigenstates of 
Td that are closest in energy to some target energy, the 
exponentially small energy gaps mean that we could be 
picking very different eigenstates (as labeled by their 1- 
bit quantum numbers) at each step. This will, in general, 
result in a slow convergence and/or a final state that is 
a superposition of many nearby eigenstates. 
Comparison with ED for small systems: We now 
benchmark our method against the Heisenberg model 
with random z-directed magnetic fields: 

H=jY,Sn-Sn+l-Y.hr,S:,. ( 2 ) 

n n 



where J = 1, are spin 1/2 operators and the fields 
are drawn randomly from the interval [—IF, IF]. This 
model has been studied extensively in the context of 
MBL and several numerical studies [7, 10, 42] indicate 
that H is fMBL for IF > 3.5 — 4. At strong disor¬ 
der, typical eigenstates look like product states in the 
CT* basis with small fluctuations. Equivalently, the “1- 
bits” t/ look like erf with exponentially decaying correc¬ 
tions from operators away from site i. As the disorder 
is lowered, the probability of many-body “Mott-type” 
resonances[43, 44], wherein the eigenstates are approxi- 
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mately equal-weight superpositions of a few basis states, 
increases. These resonant states have energy splittings 
that decay exponentially with the maximum distance in¬ 
volved in the resonance. At even smaller Ws, the ap¬ 
proach to the ergodic transition is marked by a “Grif- 
hths” region [12, 43] in which locally ergodic/critical in¬ 
clusions start to proliferate. 

Figure 1 shows a comparison between eigenenergies ob¬ 
tained using exact diagonalization (blue) and variational 
DMRG-X (red) for a system of size L = 12, disorder 
strength W = S and bond dimension y = 16. To obtain 
the full spectrum, we feed the algorithm all possible 
product states as initial states. We find that the vari¬ 
ance in the energy of all variationally obtained DMRG-X 
eigenstates is less than machine precision (~ 10“^^), and 
the overlap of these states with the exact eigenstates is 
unity up to machine precision. 

The zoomed in energy levels show that the method 
successively resolves the exponentially small splittings in 
the spectrum extremely accurately. However, a few ex¬ 
act eigenenergies have no DMRG-X partner—when two 
or more eigenstates of H have maximum weight on the 
same input basis state, the input state converges to one 
of these eigenstates leaving the other unpaired. We can 
avoid this duplication by requiring every new state to be 
orthogonal to the prior ones, but this will not be nec¬ 
essary in larger systems where our goal will never be to 
obtain every eigenstate. 

One might worry that this is method biased towards 
product states and fails to capture resonant eigenstates. 
The bottom-right panel of Fig. 1 shows a representative 
eigenstate with a many-body “Mott” resonance involving 
a few distant basis states which is exactly captured by the 
variational state. We emphasize that the algorithm only 
uses a product state as an initial input; after that, the 
algorithm converges to the previously chosen eigenstate 
of H. As long as the bond dimension y is sufficiently large 
for the eigenstates of % to capture resonances, it is easy 
for the algorithm to converge to a resonant superposition 
starting from one of the product states with significant 
weight in the resonant eigenstate. 

Larger systems: We now turn to an evaluation of the 
algorithm for system sizes inaccessible to ED by examin¬ 
ing typical variances in the energy, = (iF^) — {H)'^ for 
approximate eigenstates of the Hamiltonian (2) obtained 
using DMRG-X at different disorder strengths W and 
system sizes L. Figure 2(a) shows the disorder averaged 
value of logj^g as a function of bond-dimension y for 
randomly chosen excited states from 200-1000 disorder 
samples at different values of W and L. The grey line at 
10“^^ marks the approximate value of machine-precision 
and we average log to capture typical behavior instead 
of deviations due to rare eigenstates. 

At strong and moderate disorder {W = 8,12), we 
see an initial rapid decrease of tr^ with y followed by a 
saturation—the saturation is expected to happen when 
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FIG. 2. (a) Disorder averaged logarithm of the energy vari¬ 

ance (7^ plotted against bond-dimension y for different disor¬ 
der strengths W and system sizes L. We see a rapid decrease 
of the typical with y for moderate-strong disorder, with 
variances falling below machine precision (shaded grey region 
below 10“^^) at small y <C 2^^^. (b), (c) Variance and { 0 - 12 ) 
plotted against DMRG time-steps for a typical run in which 
y is successively increased after each 25 steps to obtain a sin¬ 
gle eigenstate with L = 24, W = 12 using our overlap-based 
DMRG-X method (red), and a more naive energy-targeting 
method (blue) showing vastly better convergence for the over¬ 
lap method. 

the bond dimension becomes large enough to capture en¬ 
tanglement over a correlation length log 2 y ~ C- Even at 
moderate disorder W = 8, the bond dimension saturates 
quickly and y ^ 40 ^ 2^/^ is already sufficient to cap¬ 
ture states to machine precision accuracy! In this regime, 
this method can be used to really push the boundaries on 
the system sizes that we have been able to study through 
ED. As the transition to the ergodic phase is approached, 
locally thermal Griffiths regions become more probable 
and the eigenstates become more entangled. We see that 
the accuracy of the method for the small bond dimen¬ 
sions considered starts to break down around W = 5, 
though a rough extrapolation suggests that we can still 
make significant improvements by using larger y. The 
increase in variance with system-size at fixed y is to be 
expected since even clean ground state DMRG methods 
make a constant error per unit length and yield a vari¬ 
ance that grows with system size. 

At even larger sizes or/and small W, inevitable locally 
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thermal/ critical Griffiths inclusions will require special 
handling as a subset of the 1-bits now look more delocal¬ 
ized and the eigenstates have a very different structure 
from product states within the inclusions. A compar¬ 
ison with ED (not shown) on a system with an artifi¬ 
cially engineered thermal inclusion shows that the vari- 
ationally obtained states correctly capture local observ¬ 
ables away from the inclusion, but make superpositions 
between eigenstates that differ primarily in the inclu¬ 
sion region. In principle, it is possible to purify these 
states to obtain an eigenstate by using a hybrid energy- 
overlap method; an inclusion of length I occurs with 
probability [43] with < 1 and has a level spacing 
A ^ 2“^. We identify the Griffiths inclusion by looking 
for a diminished value of the frozen moment |(o'/)| in the 
states obtained by DMRG-X. We then feed these states 
into a hybrid algorithm which picks states at a chosen 
energy from the subset of states which have large over¬ 
lap with the starting state away from the inclusion. This 
ensures that we’re only trying to resolve the larger level 
spacing A <C 2“^, while also maintaining the integrity 
of the state away from the inclusion. 

Note that for a typical cut somewhere along the chain, 
the entanglement entropy scales an an area law with co¬ 
efficient ^ proportional to the localization length since 
the state looks thermal on length scales shorter than ^ 
[15]. This implies that the typical bond dimension scales 
exponentially with On the other hand, the maxi¬ 
mum entanglement entropy across all cuts in the chain 
scales logarithmically with L (a thermal region of size £ 
is exponentially rare in £, but has 0{L) chances for oc- 
curing somewhere in a system of size L, thereby giving 
£ ^ log(L)). This implies a polynomial scaling with L 
for the maximum bond-dimension Xmax- Gombining this 
with the scaling of the cost of the DMRG-X algo¬ 
rithm means that the algorithm scales exponentially in ^ 
and linearly in L if the maximum x is fixed at some 0(1) 
number ^ e^. On the other hand, if the bond dimension 
is allowed to grow to achieve a certain accuracy, then the 
cost scales polynomially in L with a power larger than 1 
and dependent on W. 

We end with two comments. First, for large system 
sizes, we can randomly sample from the spectrum and ap¬ 
proximate the underlying density of states by randomly 
choosing intial product states. Even though the DMRG- 
X sweep does not use energy targeting, we can still effec¬ 
tively target different energy densities. Deep in the dis¬ 
ordered phase, the initial product state is exponentially 
close to an actual eigenstate, and thus (H) is almost con¬ 
stant during a run. 

Second, since the variationally obtained states are 
MPSs, few-point observables can be computed extremely 
efficiently. In Fig 3(b) we show {a\af\a) plotted against 
Ea = {a\H\a) for ^ 500 variationally obtained eigen¬ 
states I a) of an MBL Hamiltonian with W = 8 and 
L = 18, 24. We see a clear violation of ETH since the 
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FIG. 3. (a) (ajcrfja) plotted against Ea = {a\H\a) for an 

ergodic system with IT = 1.5 and eigenstates ja) obtained via 
ED. This is an ETH obeying phase where the observable varies 
smoothly with energy and the fluctuations decrease with in¬ 
creasing L. (b) Same quantity evaluated with ~ 500 varia¬ 
tionally obtained eigenstates ja) of an MBL Hamiltonian with 
W = 8 and L = 18, 24 showing a clear violation of ETH 

local observable does not vary smoothly with Ea and the 
fluctuations do not decrease with L. This also lends ad¬ 
ditional support that our method is correctly capturing 
MBL eigenstates since the violation of ETH would have 
been much weaker if the states |a) were superpositions 
of actual eigenstates. Such a test is especially useful at 
large Ls where the average level spacing is smaller than 
machine precision and we need to rely on methods other 
than the variance to diagnose the goodness of the vari¬ 
ational states. Figure 3(a) shows the analogous calcula¬ 
tion in an ergodic system with IT = 1.5 and eigenstates 
obtained via ED. Here we do see a smooth variation of 
the observable with Ea, and the characteristic decrease 
in fluctuations [45] with increasing L. 

Energy targeting: We now compare the convergence 
of our overlap method with the simplest energy target¬ 
ing method which picks the eigenstate of H closest to 
a chosen energy. Figure 2(b) shows a typical DMRG-X 
run with L = 24 and IT = 12 to obtain a single eigen¬ 
state. The bond dimension is increased every 25 steps 
and takes the values y = (4,8,16,24,32). We see that 
the overlap method (shown in red) converges extremely 
quickly each time the bond-dimension is increased and 
rapidly reaches machine precision. On the other hand, 
the energy targeting method (shown in blue) run for the 
same disorder realization and a target energy equal to 
the energy of the state obtained via the overlap method 
(upto 4 digits of precision) shows an extremely poor con¬ 
vergence and very large variances. In Fig. 2(c) we plot 
the expectation value of for a site in the middle of the 
chain evaluated using the states at each DMRG step. As 
expected, the overlap method shows very little fluctua¬ 
tion in this quantity, while the naive energy approach is 
clearly seen to be rattling between states with extremely 
different local quantum numbers. 

Summary and Outlook: In summary, we have devel¬ 
oped a DMRG-X method that successfully obtains highly 
excited eigenstates of MBL systems to machine precision 
accuracy at moderate-large disorder in a time that scales 
only polynomially with L. This method explicitly takes 
advantage of the local spatial structure and order char- 
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acterizing MBL eigenstates, thereby moving away from 
traditional energy based DMRG algorithms. 

A natural next step is to use the DMRG-X method to 
obtain phase boundaries between localized phases with 
different kinds of eigenstate order present. The nature of 
the phase transition between different localized phases is 
an important open question, and refining this technique 
to access these transitions at larger system sizes should 
help settle some of these questions. 
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METHOD 

In this section, we briefly recapitulate the standard DMRG algorithm[l] implemented in the language of matrix- 
product states (MPSs) [2], before providing more details on the DMRG-X method. Section IV in Ref. 3 also provides 
a very clear exposition of these numerical methods and some of our discussion closely follows this reference. 

A general quantum state j'k) for a one-dimensional system of L sites can be written in the following matrix-product 
state (MPS) form: 

1'*^)= E E (1) 

0<7„<=x,> 

where |j„) with = 1,..., d is a basis of local states at site n (for a spin 1/2 system, d = 2 and |j„) = | f), | i)), and 
the are rank three tensors (except on the first and last sites where they are rank two tensors). Figure 1(a) shows 
a pictorial representation of an MPS. The enternal legs are the “physical” spin indices whereas the internal legs 
7 „ are the virtual indices that are contracted. Each is a x„ x x„+i matrix (at the boundaries xi = xl+i = 1) 

and each matrix product ]/[^ ijWji in Eq. (1) produces a complex number which is the amplitude of |'I') on the basis 
state |ji • • • jl). 

The maximum dimension y of the {RI"!} matrices is called the bond-dimension of the MPS and low entanglement 
states can be efficiently represented my MPSs of bond dimension y <C The relationship between y and the 

entanglement can be made more precise by considering the Schmidt decomposition of the state |'I'). For a given 
bipartition of the system into left and right halves, a singular value decomposition can be used to rewrite j'k) as 

1^) = Aa|a)L|Q:}fl 
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FIG. 1. Diagrammatic representation of (a) the state |v[') as an MPS, (b) !<&) as a canonical MPS, (c) the Hamiltonian H as 
an MPO, (d) co-efficients 0 of It?) in the variational basis and (e) the effective Hamiltonian T-i in the variational basis. 
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where the form orthonormal bases for the left and right halves respectively, l(q;|/3)l = _r(q;|/?)_r = and 

the entanglement entropy of the bipartition is defined through the Schmidt values Aq as Se = — l^aPlnjAap. 
Following a prescription by Vidal[4], it is possible to define a canonical form (Fig.1(b)) for the MPS by rewriting each 
matrix as a product of a x„ x Xn+i dimensional complex matrix fI"!-!" and a square diagonal matrix A^"! such 

that matrices Al"! matrices contain the non-zero Schmidt values for a bipartition between sites n and n -I- 1 

|vE-)= ^ rW^'^AWF[2l^-^A[2l...A[^-ilF[^l^-|jd,...,ji) = |^ AN|a„)L|a„)i^, (2) 


and the states 


|a„)i= ^ (FW^-^AW..-A[”-ilr”)„|ji---j„) 

!«„)«= ^ (F["+il^"iA["+il--.A[^-ilr^)„|j„+i.-.j,.) (3) 

3n + l---jL 

define the orthonormal Schmidt states for the left and right halves of the bipartition respectively. This canonical form 
clearly relates the bond dimension of the MPS x to the number of Schmidt values contributing signihcantly to the 
entanglement entropy. 

A standard two-site DMRG algorithm tries to find the ground state by variationally optimizing the MPS 
matrices on neighboring sites to minimize the energy {ifQ\H\ipf) while keeping the rest of the chain fixed. We define 
a matrix-product operator (MPO) representation of H exactly analagous to Eq. 1 but now using 4-index tensors M 
as shown in Fig. 1(c), and then implement the following steps: 

• To update the MPS between sites n and n -I- 1, rewrite the state [T) in the basis spanned by the states |j„), 
\kn+i) and the left and right Schmidt states |a„_i)L and \ j3n+i)R 

IV') = X] ^^ap\0‘ri-l)L\jn)\kn+l)\l3n+l) 

j,k,a,(3 

where the definition of 0 is shown pictorially in Fig. 1(d) and follows directly from the definitions (3). 

• Dehne an effective Hamiltonian H as the Hamiltonian projected to the \ajkp) basis. This is a (Px^ ^ 
dimensional operator, again best seen pictorially in Fig. 1(e). 

• Find the ground state of "H, denoted by 0^^- This is the optimal state for minimizing {ifo\H\ipQ) in this subspace. 

• Do an an SVD on 0 to put the MPS back in canonical form with the matrices pM^Xl"!,updated. 

• Repeat for the next pair of sites and iteratively sweep through the chain till the state converges. 

The only difference between the ground-state DMRG algorithm outlined above and DMRG-X is in step 3. In 
the DMRG-X algorithm, we find all the Px^ eigenstates of H instead of just its ground state. We then pick 0 
as the eigenstate of H with the maximum overlap with the previously found state in the iterative scheme. The 
algorithm is initialized with an appropriate initial state which is perturbatively “close” to the true eigenstates 
of the MBL Hamiltonian (such as a product state in the cr^ basis). 
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